The goal of study 2 was to…
library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
source("functions.R")
df <- read.csv("../data/study2.csv") |>
mutate(
Date = lubridate::dmy(Date),
Participant = fct_reorder(Participant, Date),
Screen_Refresh = as.character(Screen_Refresh),
Illusion_Side = as.factor(Illusion_Side),
# Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Null", "Congruent"),
Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Congruent"),
Block = as.factor(Block),
Education = fct_relevel(Education, "Doctorate", "Master", "Bachelor", "High School", "Other")
)
# Fix precision
for(ill in unique(df$Illusion_Type)) {
data <- df[df$Illusion_Type == ill, ]
i <- 10
while (length(sort(unique(signif(data$Illusion_Difference, i)))) != 8) {
i <- i - 1
}
df[df$Illusion_Type == ill, "Illusion_Difference"] <- signif(df[df$Illusion_Type == ill, "Illusion_Difference"], i)
if (i != 10) {
message(ill, ": Illusion_Difference values != 8. Rounded to ", i, ".")
}
}
# Transformation
df <- df |>
mutate(
Illusion_Difference_log = log(1 + Illusion_Difference),
Illusion_Difference_sqrt = sqrt(Illusion_Difference),
Illusion_Difference_cbrt = Illusion_Difference**(1 / 3),
Illusion_Strength_log = sign(Illusion_Strength) * log(1 + abs(Illusion_Strength)),
Illusion_Strength_sqrt = sign(Illusion_Strength) * sqrt(abs(Illusion_Strength)),
Illusion_Strength_cbrt = sign(Illusion_Strength) * (abs(Illusion_Strength)**(1 / 3))
)
# plot(estimate_density(dfraw[dfraw$Participant == "60684f29dbfe1bb2059e5e27_rkqoy", "RT"]))
# Dear participant, thank you for participating in our studfy. Unfortunately, we didn't receive your data :( did something happen? Some technical issue? We would like to kindly ask you to return your participation so that we can open up more slots. Thank you in advance, and apologies for the inconvenience!
# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers), which makes it unusable. We understand that you might have been in a hurry or had some other issues; we hope to open-up more slots in the future would you be interested to participate again.
# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers - in particular in the 2nd part of the study), which makes it unusable. We understand that you might have been in a hurry or had some other issues, and so we kindly ask you to return your participation; we hope to open-up more slots in the future would you be interested to participate again.
# Just received the results: in your case, the three most prominent issues were your response pattern that was equivalent to random (in the sequence of keystrokes) that led to a chance-level performance (that was also significantly different from the rest of the population). Moreover, your reaction time distribution was also very different from the norm, with a vast majority of implausibly short responses (i.e., that are faster than the time it takes the brain to process a visual input). This issue was even more prominent in the second block (after the break), which typically happens when participants are in a rush to finish. Finally, your overall completion time was also significantly below the average. Again, we apologize, we understand that your time is valuable, but unfortunately we run too on limited funds and can hardly spare more on unusable data. Sorry for the inconvenience, we will try to open-up more slots soon.
outliers <- c(
# Error rate of 48.8% Very short RT
# Prolific Status: REJECTED (06/08)
"5e66ed53de7896486f9a71f8_p8ckk",
# 2nd block of responses very fast
# Prolific Status: REJECTED (15/08)
"61443e47c248a375c899a0cf_9qz3u"
)
partial_outliers <- c(
# 2nd block a bit bad
"5c43cd414fe4f800016e4983_0zp37",
# Entire 2nd block bad
"61564e378e974bdbe42763a2_hhufm")
We removed 2 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.
For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).
dfsub <- df |>
group_by(Participant) |>
summarize(
# n = n(),
Error = sum(Error) / n(),
RT_Mean = mean(RT),
RT_SD = sd(RT),
IPIP6_SD = mean(IPIP6_SD),
PID5_SD = mean(PID5_SD),
) |>
ungroup() |>
arrange(desc(Error))
data.frame(Participant = c("Total"), t(sapply(dfsub[2:ncol(dfsub)], mean, na.rm=TRUE))) |>
rbind(dfsub) |>
knitr::kable() |>
kableExtra::row_spec(1, italic = TRUE) |>
kableExtra::row_spec(which(dfsub$Participant %in% outliers) + 1, background = "#EF9A9A") |>
kableExtra::row_spec(which(dfsub$Participant %in% partial_outliers) + 1, background = "#FFECB3")
| Participant | Error | RT_Mean | RT_SD | IPIP6_SD | PID5_SD |
|---|---|---|---|---|---|
| Total | 0.220 | 793 | 601 | 20.34 | 20.48 |
| 5e66ed53de7896486f9a71f8_p8ckk | 0.488 | 311 | 337 | NA | NA |
| 61443e47c248a375c899a0cf_9qz3u | 0.359 | 417 | 434 | 19.93 | 21.55 |
| 60d60bb9ced0566454d42750_020cx | 0.349 | 618 | 227 | 22.37 | 25.55 |
| 5c43cd414fe4f800016e4983_0zp37 | 0.339 | 796 | 657 | 23.26 | 12.29 |
| 60d84da28f38998c15e2a4fa_xw2cl | 0.333 | 702 | 376 | 25.12 | 31.38 |
| 5cbb77bf418c540001641e67_sp53j | 0.324 | 870 | 547 | 23.39 | 15.09 |
| 60c5ae1d153847577f543609_u9w9l | 0.302 | 878 | 950 | 32.58 | 46.24 |
| 611d69ffd22b95ca500604b2_hs0kt | 0.300 | 855 | 528 | 24.03 | 15.30 |
| 6140b5b201e55201433c2383_60sey | 0.299 | 674 | 303 | 9.79 | 16.05 |
| 5f7cd5bc00f2a016f22e74d2_4qe89 | 0.296 | 521 | 269 | 25.43 | 36.14 |
| 61564e378e974bdbe42763a2_hhufm | 0.289 | 663 | 645 | 25.81 | 25.94 |
| 5c067e6329abc0000199ae45_lbvy6 | 0.283 | 579 | 178 | 14.20 | 17.66 |
| 60f1ce9c530c91c0a6e9dab0_1rwja | 0.283 | 766 | 399 | 22.55 | 23.73 |
| 5efbc3a84e1b6416ad27a3ec_qdax2 | 0.281 | 837 | 2198 | 14.87 | 11.78 |
| 60e72879166656b7391b2eb2_4t9nj | 0.280 | 897 | 512 | 19.71 | 25.95 |
| 6059025fd396ddcd93fe50aa_pytn9 | 0.279 | 1112 | 1363 | 20.25 | 23.23 |
| 613bb8d7741ae6ca04e0431b_gk94a | 0.277 | 921 | 565 | NA | NA |
| 5cb632bab9c4860001700fe6_dctvu | 0.277 | 719 | 570 | 16.66 | 14.08 |
| 613d8489ab93bd9635c12051_zwbj3 | 0.272 | 730 | 269 | 31.39 | 23.54 |
| 60560efe16c645d2a9bd0daa_llvby | 0.270 | 744 | 301 | 22.17 | 18.34 |
| 613c1ac80e6fc9ac43eac753_dbdal | 0.266 | 743 | 396 | 20.36 | 18.70 |
| 6170adf5330be4e266555619_68ago | 0.262 | 824 | 542 | 18.08 | 13.53 |
| 614b55e22ff3944a165736bb_v16g9 | 0.261 | 518 | 170 | 11.49 | 14.03 |
| 5eaef8702b68455d6e130595_wh7q8 | 0.259 | 707 | 338 | 14.92 | 12.77 |
| 606781fa7584b511b845d52c_cu5nv | 0.256 | 550 | 152 | NA | NA |
| 5f9aba6600cdf11f1c9b915c_at42y | 0.254 | 788 | 593 | 22.87 | 20.80 |
| 5f2e1510b83f09172c62ba16_22y38 | 0.251 | 1023 | 1708 | 8.40 | 16.45 |
| 60feee53991033aa986fa755_lp5dv | 0.251 | 658 | 330 | 13.83 | 17.35 |
| 5bdbfff9ee652a0001efca64_vvdyu | 0.248 | 627 | 302 | 23.10 | 19.36 |
| 61438c83c0d9428ad7b4e4b8_t0wjb | 0.248 | 731 | 336 | 22.27 | 34.83 |
| 612383021976c4e2cac79581_0fzdd | 0.247 | 738 | 235 | 30.22 | 23.96 |
| 60087e9189ddb34b79b08be5_pd14q | 0.245 | 613 | 202 | 26.64 | 18.85 |
| 589630d22a697d0001cfca12_buz00 | 0.243 | 850 | 253 | 18.60 | 30.58 |
| 5f46a8269f9d230534df642a_gql1v | 0.242 | 1354 | 1254 | 12.23 | 10.61 |
| 5c437f6a4fe4f800016e3d52_a73wb | 0.241 | 554 | 177 | 30.65 | 25.55 |
| 611fb1635c5f446a6bb2848d_6fcvf | 0.241 | 990 | 513 | 37.35 | 28.90 |
| 5f494676e0837811e94a4007_qfzu1 | 0.240 | 940 | 670 | 16.18 | 21.86 |
| 5f983f3b00bbab0c4b28149a_l7hnq | 0.240 | 646 | 225 | 11.55 | 21.11 |
| 5d3c6e745602310001bca8aa_6bdtb | 0.238 | 662 | 495 | NA | NA |
| 61645edbf8a9840feeb735b6_qot1q | 0.238 | 1410 | 1660 | 17.48 | 16.72 |
| 614060a52d7c64c27ef9887c_43ppx | 0.238 | 794 | 628 | 18.33 | 13.04 |
| 611e6b43e936390a65656ecd_rr1n5 | 0.238 | 578 | 166 | 27.05 | 19.75 |
| 5b8c8e7d0c740e00019d55c3_ejsw9 | 0.236 | 674 | 313 | 15.44 | 17.20 |
| 5e98b398e4a84f0eb3755e51_yy8z9 | 0.235 | 574 | 241 | 13.95 | 23.42 |
| 5c8fcb13dd3b360015b4dff6_dvbq9 | 0.234 | 643 | 327 | NA | NA |
| 5f6b027e6eae971e2fa13594_5qh7m | 0.233 | 555 | 180 | 16.80 | 13.21 |
| 5f304d07956a3f48639e3aa6_03qfs | 0.231 | 698 | 286 | 26.22 | 27.93 |
| 615d5311a9fbf1f98aa62a5a_n2x53 | 0.229 | 495 | 182 | 18.67 | 20.43 |
| 610cda1332fb63830158c55d_95zqk | 0.228 | 691 | 458 | 24.55 | 24.63 |
| 60ec51c51a3158a50ded8a3e_jj888 | 0.227 | 676 | 252 | 14.28 | 10.15 |
| 613bab4160cd2fe41f9b1a27_kepjv | 0.226 | 587 | 292 | 20.87 | 22.09 |
| 60efc61b7b9f908aea94ed57_z9d93 | 0.224 | 1506 | 12083 | 26.40 | 17.62 |
| 612cb8236aa1cef2599e6f59_6jam2 | 0.221 | 622 | 751 | 18.93 | 16.50 |
| 60d1224bb86fedf1819ef67b_tvh43 | 0.221 | 790 | 579 | 17.25 | 26.71 |
| 6166c19fb69f46ba9749489b_6ym17 | 0.220 | 735 | 313 | 22.95 | 25.24 |
| 60142d79e4c5e23b4730a69e_98qsp | 0.220 | 618 | 189 | 17.14 | 19.61 |
| 606f598baf5e0ffda45e9e47_qayzo | 0.220 | 1327 | 1557 | 9.45 | 14.82 |
| 611ee495b12c10171def6816_9oz56 | 0.219 | 1226 | 928 | 26.97 | 25.54 |
| 612819d5c4f57a5fb23d03b3_sc9zu | 0.218 | 630 | 467 | 14.06 | 16.49 |
| 611183c67911657721712872_sdh4z | 0.218 | 846 | 529 | 14.36 | 12.28 |
| 60fa8b48cbeec7cf62317c21_glept | 0.216 | 718 | 346 | 27.40 | 36.23 |
| 5c1bb460a05a64000125c522_ohxap | 0.216 | 564 | 169 | 30.73 | 20.37 |
| 60feabd18109e08e594540f8_r9m9e | 0.216 | 622 | 191 | 23.11 | 29.01 |
| 5e8e06cbdf10df0a811b6f29_g2tm6 | 0.215 | 610 | 293 | NA | NA |
| 6112810f76d46087c3e2d208_51bwf | 0.215 | 770 | 663 | 22.42 | 20.99 |
| 5e283bdff5bb81010fcc6cde_9zgbc | 0.214 | 1446 | 1197 | 13.93 | 22.34 |
| 5d2508273854c300016c04aa_1qm3y | 0.214 | 757 | 353 | 11.82 | 13.46 |
| 5b232f6838fc0c000131438c_g758c | 0.212 | 1024 | 766 | 28.16 | 21.16 |
| 5e938e072a4ce9077f8157b9_egf3b | 0.212 | 1012 | 781 | 16.75 | 27.39 |
| 5c42329384ef0f00011882ca_tae13 | 0.209 | 633 | 268 | 27.45 | 27.36 |
| 60c0fe89173c2b1f4fe2c578_a23gw | 0.209 | 783 | 341 | 31.55 | 30.88 |
| 60f5fa488f7b1381d175ecb5_qwvfs | 0.208 | 670 | 341 | 20.86 | 14.25 |
| 6010b797bbe6be440425665c_8qarc | 0.207 | 816 | 521 | 11.95 | 15.38 |
| 605d5ffca1324379006a3599_703q5 | 0.207 | 512 | 139 | 11.48 | 11.99 |
| 5e6e14f53c76d23b934a67f3_ghxwv | 0.206 | 911 | 2764 | 33.56 | 27.03 |
| 6139133041546096188da709_kr98d | 0.206 | 799 | 465 | 21.71 | 19.13 |
| 5fcbd87e70a25c99595f012d_dkj0n | 0.205 | 564 | 155 | 17.76 | 20.20 |
| 61602320c1aec7e6f8b5e3f4_adnrd | 0.205 | 558 | 177 | 22.12 | 8.45 |
| 60c21f5f1389f65d2d88a3bf_lxzqo | 0.204 | 1030 | 466 | 35.66 | 20.25 |
| 5f264332488e13372281c1c8_b082g | 0.204 | 745 | 359 | 19.00 | 11.51 |
| 60a5404314ae1a9d5b88565d_fd4x8 | 0.204 | 763 | 274 | 15.46 | 15.63 |
| 61264d8a5de54eb1e42ed381_6x9p2 | 0.203 | 735 | 242 | 25.56 | 6.79 |
| 60009275300e1c14bf67ce83_7n71t | 0.203 | 602 | 202 | 22.41 | 31.03 |
| 5fa18c4a48882016b4b143f7_kn4c0 | 0.202 | 965 | 984 | 20.97 | 32.36 |
| 62d6d791bb6448ae52929ebd_ncqkx | 0.202 | 776 | 467 | 30.03 | 27.75 |
| 61395b705a38ede105ed7982_jkfrp | 0.202 | 639 | 298 | 20.92 | 25.94 |
| 5f439ac86b921f5fac4f55ad_afpfo | 0.196 | 884 | 656 | 22.10 | 16.21 |
| 60127f78468af02bcca0b7a0_lad8d | 0.196 | 576 | 250 | NA | NA |
| 5f97e6601f6d0e016087fc91_bac1g | 0.195 | 649 | 223 | 20.92 | 19.93 |
| 614f8321333460db43b7f229_mpxhl | 0.194 | 1645 | 983 | NA | NA |
| 5f733312e59fd2000a74b8cf_lbgw8 | 0.194 | 972 | 612 | 12.46 | 16.09 |
| 5fc78656c7368206ab1c5174_4t6uk | 0.192 | 640 | 306 | 13.37 | 22.43 |
| 614b23f447bd48a3801fbf95_dlqf9 | 0.190 | 739 | 437 | 16.14 | 18.70 |
| 6106ac34408681f3b0d07396_y2n16 | 0.189 | 613 | 221 | 27.84 | 24.74 |
| 60e866921fd3ba5f3007c2aa_7toga | 0.189 | 987 | 1084 | 16.00 | 23.38 |
| 5ed7a7a467a98224295459ff_0l6hx | 0.187 | 531 | 107 | 13.22 | 10.44 |
| 5ceef0cd45c1ec001920d548_kurpx | 0.185 | 1062 | 758 | 22.10 | 22.66 |
| 610af9910c4d079d84e855f4_z8j00 | 0.185 | 764 | 360 | 25.17 | 16.11 |
| 5e6d381c029fa42f9f48c105_kcgfo | 0.185 | 869 | 562 | 24.09 | 20.61 |
| 5f782504951dda38db3ef0b6_qsfgg | 0.184 | 1072 | 798 | 14.89 | 15.58 |
| 613646a48d30fb040cf5ffed_g5q2e | 0.184 | 702 | 184 | 14.16 | 15.04 |
| 614a4c9d8172751921452626_mc4l4 | 0.184 | 772 | 305 | 31.56 | 28.66 |
| 5fdfd04b9bf07d83b2e5f780_qcv0k | 0.183 | 886 | 537 | 23.91 | 23.77 |
| 613231ee733769a9aff9c6c5_ps1tt | 0.180 | 896 | 619 | 11.43 | 25.06 |
| 60d26bfdacc7efc851571e74_l8k5r | 0.180 | 718 | 372 | 18.15 | 27.52 |
| 5eb3a734d249ac18a413063a_ec0gy | 0.179 | 921 | 527 | 20.16 | 14.07 |
| 5f61341fe5033f0b6327fb7f_lg9vw | 0.179 | 775 | 384 | 13.89 | 20.64 |
| 5ede636d21cf560bf7aeb29f_03zzd | 0.178 | 741 | 532 | 17.33 | 21.62 |
| 6026f1c19366fc4778bfd055_nuccw | 0.173 | 873 | 568 | 15.28 | 19.65 |
| 60ddfb3db6a71ad9ba75e387_9elmo | 0.173 | 692 | 287 | 16.89 | 12.46 |
| 60d340fe43207d6a5bc89ca1_m32vm | 0.173 | 1090 | 1057 | 14.99 | 18.71 |
| 5f31de301a469e11624f7948_3n2lv | 0.164 | 783 | 311 | 19.53 | 26.15 |
| 611bb2a6d9678117db737236_pf2qy | 0.163 | 904 | 549 | 12.64 | 12.23 |
| 6147c21c933739516362d0ee_wcaat | 0.163 | 900 | 460 | 24.99 | 24.26 |
| 5ecd58b69c382e0a177e0f61_p3yev | 0.162 | 652 | 279 | NA | NA |
| 6067f370472020b9dc5cf71b_f2or0 | 0.162 | 861 | 337 | 16.49 | 17.12 |
| 61034f24da19cc56177b8b59_43rrv | 0.161 | 805 | 418 | 28.19 | 13.51 |
| 614729093372b9e84bc48f8e_v4mcv | 0.159 | 952 | 530 | 26.37 | 17.89 |
| 62da72f3fb2ff174670bc85f_vb5a9 | 0.157 | 786 | 365 | 18.68 | 13.44 |
| 5e8348450e128b0539d77ddc_z361c | 0.154 | 820 | 523 | 9.39 | 24.93 |
| 5f8cbc5c355ea745e6cef2ca_5wo4f | 0.152 | 765 | 406 | 10.78 | 23.51 |
| 610a9a0723dac03a63f15035_gclx4 | 0.151 | 795 | 312 | 28.52 | 35.59 |
| 612684c6d2cef3acadf5d77b_pv444 | 0.151 | 668 | 389 | 26.37 | 26.56 |
| 5d6e3252dfe6d10001dbe508_1aloo | 0.150 | 816 | 379 | 25.32 | 19.75 |
| 60cefa69352cbf2549f2bf35_e13a8 | 0.150 | 688 | 221 | 24.42 | 25.38 |
| 5efd2964d36f63162f263795_xksu5 | 0.148 | 849 | 646 | 9.73 | 14.50 |
| 60ec6cb208c74c0d7e4794b9_hjay3 | 0.147 | 865 | 720 | 18.17 | 17.81 |
| 56cb8858edf8da000b6df354_m4821 | 0.146 | 579 | 211 | 12.57 | 10.43 |
| 6165e0b51a883e9db8cc7146_z0xax | 0.145 | 816 | 647 | 23.25 | 14.69 |
| 5cdc88cc50f50f001a783112_tw4fr | 0.118 | 1933 | 1570 | 23.15 | 13.20 |
# RT distribution
p <- df |>
# filter(Date > "2022-08-12") |>
filter(RT < 10000) |>
estimate_density(select = "RT", at = c("Participant", "Block")) |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
mutate(Participant = factor(Participant, levels=levels(df$Participant))) |>
mutate(color = ifelse(Participant %in% outliers, "red", ifelse(Participant %in% partial_outliers, "orange", "blue"))) |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
geom_vline(xintercept = 125, linetype = "dashed", color = "red") +
scale_color_manual(values = c("red" = "#F44336", "orange" = "#FF9800", "blue" = "blue"), guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(xlim = c(0, 3000)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
facet_wrap(~Participant) +
labs(y = "", x = "Reaction Time (ms)")
# p
# ggsave("figures/outliers.png", p, width=20, height=15)
# Filter out
df <- filter(df, !Participant %in% outliers)
temp <- df |>
group_by(Participant, Illusion_Type, Block) |>
summarize(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
arrange(desc(ErrorRate_per_block))
temp2 <- temp |>
filter(ErrorRate_per_block >= 0.5) |>
group_by(Illusion_Type, Block) |>
summarize(n = n()) |>
arrange(desc(n), Illusion_Type) |>
ungroup() |>
mutate(
n_trials = cumsum(n * 56),
p_trials = n_trials / nrow(df)
)
# knitr::kable(temp2)
p1 <- temp |>
estimate_density(at = c("Illusion_Type", "Block"), method="KernSmooth") |>
ggplot(aes(x = x, y = y)) +
geom_line(aes(color = Illusion_Type, linetype = Block)) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(y = "Distribution", x = "Error Rate") +
theme_modern()
p2 <- temp2 |>
mutate(Block = fct_rev(Block)) |>
ggplot(aes(x = Illusion_Type, y = p_trials)) +
geom_bar(stat = "identity", aes(fill = Block)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
theme_modern() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1 | p2
# Drop
df <- df |>
group_by(Participant, Illusion_Type, Block) |>
mutate(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
filter(ErrorRate_per_block < 0.5) |>
select(-ErrorRate_per_block)
# Drop also participant with bad second block
df <- filter(
df,
!(Participant %in% partial_outliers & df$Block == 2))
rm(temp, temp2)
df <- df |>
group_by(Participant, Error) |>
mutate(Outlier = ifelse(Error == 0 & (RT < 125 | standardize(RT) > 4), TRUE, FALSE)) |>
ungroup()
p1 <- estimate_density(df, select = "RT", at = "Participant") |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
merge(df |>
filter(Error == 0) |>
group_by(Participant) |>
summarize(Outlier = mean(RT) + 4 * sd(RT))) |>
mutate(Outlier = ifelse(x >= Outlier, TRUE, FALSE)) |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = Participant, linetype = Outlier), alpha=0.2) +
geom_vline(xintercept = c(125), linetype = "dashed", color = "red") +
scale_color_material_d("rainbow", guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
guides(linetype = "none") +
coord_cartesian(xlim = c(0, 4000)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
labs(y = "", x = "Reaction Time (ms)")
p2 <- df |>
group_by(Participant) |>
summarize(Outlier = sum(Outlier) / n()) |>
mutate(Participant = fct_reorder(Participant, Outlier)) |>
ggplot(aes(x = Participant, y = Outlier)) +
geom_bar(stat = "identity", aes(fill = Participant)) +
scale_fill_material_d("rainbow", guide = "none") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
see::theme_modern() +
theme(axis.text.x = element_blank()) +
labs(y = "Percentage of outlier trials")
p1 / p2
We removed 1244 (0.78%) outlier trials (125 ms < RT < 4 SD above mean).
df <- filter(df, Outlier == FALSE)
df$RT <- df$RT / 1000 # Convert to second for better model convergence
dfsub <- df |>
group_by(Participant) |>
select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS, starts_with("IPIP6"), starts_with("PID5")) |>
slice(1) |>
ungroup()
The final sample included 128 participants (Mean age = 26.1, SD = 7.1, range: [18, 69]; Sex: 46.9% females, 52.3% males, 0.8% other; Education: Doctorate, 0.78%; Master, 12.50%; Bachelor, 46.09%; High School, 34.38%; Other, 5.47%; Prefer not to Say, 0.78%).
plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
dfsub |>
ggplot(aes_string(x = what)) +
geom_density(fill = fill) +
geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
ggtitle(title, subtitle = subtitle) +
theme_modern() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(face = "italic", hjust = 0.5),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank()
)
}
plot_waffle <- function(dfsub, what = "Nationality", title = what) {
ggwaffle::waffle_iron(dfsub, what) |>
# mutate(label = emojifont::fontawesome('fa-twitter')) |>
ggplot(aes(x, y, fill = group)) +
ggwaffle::geom_waffle() +
# geom_point() +
# geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
coord_equal() +
ggtitle(title) +
labs(fill = "") +
theme_void() +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")
p4 <- plot_waffle(dfsub, "Sex") +
scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))
p5 <- plot_waffle(dfsub, "Education") +
scale_fill_viridis_d()
p6 <- plot_waffle(dfsub, "Nationality") +
scale_fill_manual(values = c("South American" = "#FF5722", "Middle Eastern" = "#FF9800", "Western" = "#2196F3", "African" = "#4CAF50", "South African" = "#009688"))
p7 <- plot_waffle(dfsub, "Ethnicity") +
scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0", "Mixed" = "#795548"))
p8 <- plot_waffle(dfsub, "Screen_Resolution", title = "Screen Resolution") +
scale_fill_pizza_d() +
guides(fill = "none")
p9 <- plot_waffle(dfsub, "Device_OS", title = "Device OS") +
scale_fill_bluebrown_d()
# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
# scale_fill_viridis_d()
(p1 / p2 / p3) | (p4 / p5 / p9) | (p6 / p7 / p8)
data <- filter(df, Illusion_Type == "Delboeuf")
# mm <- model.matrix(Error ~ abs(Illusion_Strength) * Illusion_Effect, data=data)
# head(as.data.frame(mm))
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_sqrt 11023 0.4205 6.62e-05 Errors
## 2 DIFF_cbrt--STRENGTH_sqrt 11026 0.4125 = 0.195 6.63e-05 Errors
## 3 DIFF_log--STRENGTH_sqrt 11038 0.4335 < 0.001 6.73e-05 Errors
## 4 DIFF_sqrt--STRENGTH_cbrt 11046 0.4345 < 0.001 6.91e-05 Errors
## 5 DIFF_sqrt--STRENGTH_log 11048 0.4116 < 0.001 6.62e-05 Errors
## 6 DIFF_sqrt--STRENGTH_log 11418 0.0710 2.70e-01 RT
## 7 DIFF_log--STRENGTH_log 11421 0.0709 = 0.175 2.66e-01 RT
## 8 DIFF_cbrt--STRENGTH_log 11422 0.0708 = 0.122 2.72e-01 RT
## 9 DIFF_sqrt--STRENGTH_sqrt 11432 0.0703 < 0.001 2.60e-01 RT
## 10 DIFF_sqrt--STRENGTH 11433 0.0702 < 0.001 2.72e-01 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_delboeuf_err <- brms::brm(formula,
data = data,
refresh = 0
)
# parameters::parameters(model_delboeuf_err)
save(model_delboeuf_err, file="models/model_delboeuf_err.Rdata")
load("models/model_delboeuf_err.RData")
p_delboeuf_err <- plot_model_err(data, model_delboeuf_err)
p_delboeuf_err
data <- filter(df, Illusion_Type == "Delboeuf", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
sigma ~ 1 + (1 | Participant),
beta ~ 1 + (1 | Participant),
family = "exgaussian"
)
# brms::get_prior(formula, data = data)
# plot(seq(-5, 5, length.out=100), brms::dstudent_t(seq(-5, 5, length.out=100), df = 3, 0.1, 0.5))
priors <- c(
# Intercept
prior("student_t(3, 0, 0.1)", class = "Intercept"),
prior("student_t(3, 0.35, 0.1)", class = "sd", coef="Intercept", group="Participant"),
# Default parameters
prior("student_t(3, 0, 3)", class = "b"),
prior("student_t(3, 0.35, 0.5)", class = "sd", group="Participant")
# prior("student_t(3, 0, 1)", class="Intercept", dpar = "sigma"),
# prior("student_t(3, 0.1, 1)", class="sd", dpar = "sigma")
# prior("student_t(3, 0.1, 0.5)", class = "ndt")
) |>
brms::validate_prior(formula, data = data)
# bayestestR::model_to_priors(model_delboeuf_rt, scale_multiply=3)
model_delboeuf_rt <- brms::brm(formula,
data = data,
refresh = 50,
init=0
)
save(model_delboeuf_rt, file="models/model_delboeuf_rt.Rdata")
load("models/model_delboeuf_rt.RData")
plot_ppcheck(model_delboeuf_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
p_delboeuf_rt <- plot_model_rt(data, model_delboeuf_rt)
p_delboeuf_rt
p_delboeuf <- plot_all(data, p_delboeuf_err, p_delboeuf_rt)
p_delboeuf
data <- filter(df, Illusion_Type == "Ebbinghaus")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_log 8184 0.5750 8.73e-09 Errors
## 2 DIFF_log--STRENGTH_log 8185 0.5948 = 0.454 9.21e-09 Errors
## 3 DIFF_sqrt--STRENGTH 8193 0.5638 = 0.010 8.77e-09 Errors
## 4 DIFF_log--STRENGTH 8198 0.5790 = 0.001 9.33e-09 Errors
## 5 DIFF_cbrt--STRENGTH_log 8200 0.5655 < 0.001 8.83e-09 Errors
## 6 DIFF_sqrt--STRENGTH 8336 0.0943 2.31e-01 RT
## 7 DIFF_sqrt--STRENGTH_log 8338 0.0942 = 0.447 2.33e-01 RT
## 8 DIFF_cbrt--STRENGTH 8339 0.0942 = 0.320 2.32e-01 RT
## 9 DIFF_cbrt--STRENGTH_log 8340 0.0941 = 0.154 2.34e-01 RT
## 10 DIFF_log--STRENGTH 8345 0.0938 = 0.013 2.32e-01 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_ebbinghaus_err <- brms::brm(formula,
data = data,
refresh = 0
)
# parameters::parameters(model_delboeuf_err)
save(model_ebbinghaus_err, file="models/model_ebbinghaus_err.Rdata")
load("models/model_ebbinghaus_err.RData")
p_ebbinghaus_err <- plot_model_err(data, model_ebbinghaus_err)
p_ebbinghaus_err
data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
sigma ~ 1 + (1 | Participant),
beta ~ 1 + (1 | Participant),
family = "exgaussian"
)
model_ebbinghaus_rt <- brms::brm(formula,
data = data,
refresh = 0,
init = 0
)
save(model_ebbinghaus_rt, file="models/model_ebbinghaus_rt.Rdata")
load("models/model_ebbinghaus_rt.RData")
plot_ppcheck(model_ebbinghaus_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
p_ebbinghaus_rt <- plot_model_rt(data, model_ebbinghaus_rt)
p_ebbinghaus_rt
data <- filter(df, Illusion_Type == "Rod-Frame")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_log--STRENGTH_log 10354 0.4875 0.000146 Errors
## 2 DIFF_sqrt--STRENGTH_log 10357 0.4893 = 0.208 0.000148 Errors
## 3 DIFF_log--STRENGTH_sqrt 10367 0.4845 = 0.002 0.000146 Errors
## 4 DIFF_sqrt--STRENGTH_sqrt 10370 0.4865 < 0.001 0.000147 Errors
## 5 DIFF_log--STRENGTH_cbrt 10374 0.4905 < 0.001 0.000149 Errors
## 6 DIFF_log--STRENGTH_log 8688 0.0900 0.076184 RT
## 7 DIFF_log--STRENGTH_sqrt 8689 0.0899 = 0.435 0.076600 RT
## 8 DIFF_log--STRENGTH 8695 0.0895 = 0.028 0.079896 RT
## 9 DIFF_log--STRENGTH_cbrt 8700 0.0891 = 0.002 0.077029 RT
## 10 DIFF_sqrt--STRENGTH_log 8701 0.0891 = 0.001 0.077129 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
(1 + logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_rodframe_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 3 finished in 1585.6 seconds.
## Chain 1 finished in 1610.1 seconds.
## Chain 2 finished in 1616.1 seconds.
## Chain 4 finished in 1618.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1607.6 seconds.
## Total execution time: 1618.8 seconds.
# parameters::parameters(model_delboeuf_err)
save(model_rodframe_err, file="models/model_rodframe_err.Rdata")
load("models/model_rodframe_err.RData")
p_rodframe_err <- plot_model_err(data, model_rodframe_err)
p_rodframe_err
data <- filter(df, Illusion_Type == "Rod-Frame", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
sigma ~ 1 + (1 | Participant),
beta ~ 1 + (1 | Participant),
family = "exgaussian"
)
model_rodframe_rt <- brms::brm(formula,
data = data,
refresh = 0,
init = 0
)
save(model_rodframe_rt, file="models/model_rodframe_rt.Rdata")
load("models/model_rodframe_rt.RData")
plot_ppcheck(model_rodframe_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
p_rodframe_rt <- plot_model_rt(data, model_rodframe_rt)
p_rodframe_rt
data <- filter(df, Illusion_Type == "Vertical-Horizontal")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_sqrt 10047 0.579 0.128 Errors
## 2 DIFF_cbrt--STRENGTH_sqrt 10054 0.578 = 0.027 0.129 Errors
## 3 DIFF_log--STRENGTH_sqrt 10061 0.581 < 0.001 0.129 Errors
## 4 DIFF--STRENGTH_sqrt 10071 0.582 < 0.001 0.129 Errors
## 5 DIFF_sqrt--STRENGTH_cbrt 10143 0.586 < 0.001 0.129 Errors
## 6 DIFF_log--STRENGTH_sqrt 7234 0.128 0.875 RT
## 7 DIFF_sqrt--STRENGTH_sqrt 7234 0.128 = 0.956 0.888 RT
## 8 DIFF--STRENGTH_sqrt 7237 0.128 = 0.189 0.871 RT
## 9 DIFF_cbrt--STRENGTH_sqrt 7238 0.128 = 0.089 0.893 RT
## 10 DIFF_log--STRENGTH_cbrt 7247 0.127 = 0.001 0.880 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_verticalhorizontal_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 1850.4 seconds.
## Chain 4 finished in 1869.3 seconds.
## Chain 2 finished in 1916.9 seconds.
## Chain 3 finished in 1929.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1891.5 seconds.
## Total execution time: 1929.7 seconds.
# parameters::parameters(model_verticalhorizontal_err)
save(model_verticalhorizontal_err, file="models/model_verticalhorizontal_err.Rdata")
load("models/model_verticalhorizontal_err.RData")
p_verticalhorizontal_err <- plot_model_err(data, model_verticalhorizontal_err)
p_verticalhorizontal_err
data <- filter(df, Illusion_Type == "Vertical-Horizontal", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
# sigma ~ t2(Illusion_Difference * Illusion_Strength) +
# (1 | Participant),
family = "exgaussian"
)
model_verticalhorizontal_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_verticalhorizontal_rt, file="models/model_verticalhorizontal_rt.Rdata")
load("models/model_verticalhorizontal_rt.RData")
plot_ppcheck(model_verticalhorizontal_rt)
# performance::check_model(model_verticalhorizontal_rt)
# parameters::parameters(model)
p_verticalhorizontal_rt <- plot_model_rt(data, model_verticalhorizontal_rt)
p_verticalhorizontal_rt
data <- filter(df, Illusion_Type == "Zöllner")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_cbrt--STRENGTH 9367 0.4069 0.71416 Errors
## 2 DIFF_log--STRENGTH 9375 0.4154 = 0.013 0.71320 Errors
## 3 DIFF_sqrt--STRENGTH 9376 0.4170 = 0.009 0.71262 Errors
## 4 DIFF--STRENGTH 9521 0.4468 < 0.001 0.70663 Errors
## 5 DIFF_cbrt--STRENGTH_sqrt 9772 0.3878 < 0.001 0.70402 Errors
## 6 DIFF_cbrt--STRENGTH 13383 0.0690 0.00455 RT
## 7 DIFF_log--STRENGTH 13384 0.0689 = 0.545 0.00467 RT
## 8 DIFF_sqrt--STRENGTH 13396 0.0683 = 0.001 0.00483 RT
## 9 DIFF_cbrt--STRENGTH_sqrt 13420 0.0670 < 0.001 0.00442 RT
## 10 DIFF_log--STRENGTH_sqrt 13420 0.0669 < 0.001 0.00455 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
(1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_zollner_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 3 finished in 1926.9 seconds.
## Chain 1 finished in 1964.4 seconds.
## Chain 2 finished in 1973.6 seconds.
## Chain 4 finished in 1980.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1961.3 seconds.
## Total execution time: 1980.3 seconds.
# parameters::parameters(model_zollner_err)
save(model_zollner_err, file="models/model_zollner_err.Rdata")
load("models/model_zollner_err.RData")
p_zollner_err <- plot_model_err(data, model_zollner_err)
p_zollner_err
data <- filter(df, Illusion_Type == "Zöllner", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
(1 | Participant),
# sigma ~ t2(Illusion_Difference * Illusion_Strength) +
# (1 | Participant),
family = "exgaussian"
)
model_zollner_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_zollner_rt, file="models/model_zollner_rt.Rdata")
load("models/model_zollner_rt.RData")
plot_ppcheck(model_zollner_rt)
# performance::check_model(model_zollner_rt)
# parameters::parameters(model)
p_zollner_rt <- plot_model_rt(data, model_zollner_rt)
p_zollner_rt
data <- filter(df, Illusion_Type == "White")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_sqrt 5821 0.7677 0.00175 Errors
## 2 DIFF--STRENGTH_sqrt 5827 0.7809 = 0.045 0.00165 Errors
## 3 DIFF_cbrt--STRENGTH_sqrt 5834 0.7644 = 0.002 0.00182 Errors
## 4 DIFF_log--STRENGTH_sqrt 5861 0.7615 < 0.001 0.00193 Errors
## 5 DIFF_sqrt--STRENGTH_log 5862 0.7970 < 0.001 0.00200 Errors
## 6 DIFF_log--STRENGTH 9590 0.0903 0.01015 RT
## 7 DIFF_cbrt--STRENGTH 9600 0.0897 = 0.008 0.01027 RT
## 8 DIFF_sqrt--STRENGTH 9610 0.0892 < 0.001 0.01035 RT
## 9 DIFF_log--STRENGTH_sqrt 9644 0.0866 < 0.001 0.01129 RT
## 10 DIFF--STRENGTH 9648 0.0869 < 0.001 0.01048 RT
plot_descriptive_err(data)
formula <- brms::bf(
Error ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_white_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Start sampling
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 2677.9 seconds.
## Chain 4 finished in 2991.8 seconds.
## Chain 3 finished in 3058.0 seconds.
## Chain 2 finished in 3117.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 2961.3 seconds.
## Total execution time: 3117.8 seconds.
# parameters::parameters(model_white_err)
save(model_white_err, file="models/model_white_err.Rdata")
load("models/model_white_err.RData")
p_white_err <- plot_model_err(data, model_white_err)
p_white_err
data <- filter(df, Illusion_Type == "White", Error == 0)
plot_descriptive_rt(data)
formula <- brms::bf(
RT ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
# sigma ~ t2(Illusion_Difference * Illusion_Strength) +
# (1 | Participant),
family = "exgaussian"
)
model_white_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_white_rt, file="models/model_white_rt.Rdata")
load("models/model_white_rt.RData")
plot_ppcheck(model_white_rt)
# performance::check_model(model_white_rt)
# parameters::parameters(model)
p_white_rt <- plot_model_rt(data, model_white_rt)
p_white_rt
data <- filter(df, Illusion_Type == "Müller-Lyer")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_sqrt 7273 0.732 5.25e-17 Errors
## 2 DIFF_log--STRENGTH_sqrt 7273 0.730 = 0.738 4.28e-17 Errors
## 3 DIFF--STRENGTH_sqrt 7291 0.730 < 0.001 4.32e-17 Errors
## 4 DIFF_cbrt--STRENGTH_sqrt 7292 0.734 < 0.001 6.54e-17 Errors
## 5 DIFF_log--STRENGTH_log 7299 0.752 < 0.001 5.76e-17 Errors
## 6 DIFF_sqrt--STRENGTH 9593 0.154 1.64e-09 RT
## 7 DIFF_cbrt--STRENGTH 9593 0.154 = 0.953 1.62e-09 RT
## 8 DIFF_cbrt--STRENGTH_sqrt 9602 0.153 = 0.007 6.02e-10 RT
## 9 DIFF_sqrt--STRENGTH_sqrt 9604 0.153 = 0.004 6.19e-10 RT
## 10 DIFF_log--STRENGTH 9605 0.153 = 0.002 1.86e-09 RT
plot_descriptive_err(data, side = "updown")
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_mullerlyer_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 3 finished in 2054.6 seconds.
## Chain 1 finished in 2112.0 seconds.
## Chain 4 finished in 2506.1 seconds.
## Chain 2 finished in 2610.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 2320.9 seconds.
## Total execution time: 2611.3 seconds.
# parameters::parameters(model_white_err)
save(model_mullerlyer_err, file="models/model_mullerlyer_err.Rdata")
load("models/model_mullerlyer_err.RData")
p_mullerlyer_err <- plot_model_err(data, model_mullerlyer_err)
p_mullerlyer_err
data <- filter(df, Illusion_Type == "Müller-Lyer", Error == 0)
plot_descriptive_rt(data, side = "updown")
formula <- brms::bf(
RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
# sigma ~ t2(Illusion_Difference * Illusion_Strength) +
# (1 | Participant),
family = "exgaussian"
)
model_mullerlyer_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_mullerlyer_rt, file="models/model_mullerlyer_rt.Rdata")
# parameters::parameters(model_white_err)
load("models/model_mullerlyer_rt.RData")
plot_ppcheck(model_mullerlyer_rt)
# performance::check_model(model_mullerlyer_rt)
# parameters::parameters(model)
p_mullerlyer_rt <- plot_model_rt(data, model_mullerlyer_rt)
p_mullerlyer_rt
data <- filter(df, Illusion_Type == "Ponzo")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_cbrt--STRENGTH 8250 0.563 0.735137 Errors
## 2 DIFF_sqrt--STRENGTH 8254 0.570 = 0.135 0.737259 Errors
## 3 DIFF_log--STRENGTH 8304 0.590 < 0.001 0.742920 Errors
## 4 DIFF--STRENGTH 8339 0.599 < 0.001 0.745206 Errors
## 5 DIFF_cbrt--STRENGTH_sqrt 8368 0.586 < 0.001 0.761796 Errors
## 6 DIFF_cbrt--STRENGTH 10267 0.118 0.000664 RT
## 7 DIFF_sqrt--STRENGTH 10279 0.117 = 0.003 0.000631 RT
## 8 DIFF_log--STRENGTH 10325 0.115 < 0.001 0.000579 RT
## 9 DIFF--STRENGTH 10357 0.114 < 0.001 0.000566 RT
## 10 DIFF_cbrt--STRENGTH_sqrt 10362 0.113 < 0.001 0.001393 RT
plot_descriptive_err(data, side = "updown")
formula <- brms::bf(
Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
(1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_ponzo_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 2 finished in 2044.8 seconds.
## Chain 4 finished in 2077.4 seconds.
## Chain 3 finished in 2093.1 seconds.
## Chain 1 finished in 2096.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 2078.1 seconds.
## Total execution time: 2097.2 seconds.
# parameters::parameters(model_ponzo_err)
save(model_ponzo_err, file="models/model_ponzo_err.Rdata")
load("models/model_ponzo_err.RData")
p_ponzo_err <- plot_model_err(data, model_ponzo_err)
p_ponzo_err
data <- filter(df, Illusion_Type == "Ponzo", Error == 0)
plot_descriptive_rt(data, side = "updown")
formula <- brms::bf(
RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
(1 | Participant),
# sigma ~ t2(Illusion_Difference * Illusion_Strength) +
# (1 | Participant),
family = "exgaussian"
)
model_ponzo_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_ponzo_rt, file="models/model_ponzo_rt.Rdata")
# parameters::parameters(model_white_err)
load("models/model_ponzo_rt.RData")
plot_ppcheck(model_ponzo_rt)
# performance::check_model(model_ponzo_rt)
# parameters::parameters(model)
p_ponzo_rt <- plot_model_rt(data, model_ponzo_rt)
p_ponzo_rt
data <- filter(df, Illusion_Type == "Poggendorff")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_cbrt--STRENGTH_sqrt 8259 0.5190 0.6545 Errors
## 2 DIFF_sqrt--STRENGTH_sqrt 8281 0.5193 < 0.001 0.6607 Errors
## 3 DIFF_cbrt--STRENGTH 8325 0.4986 < 0.001 0.6401 Errors
## 4 DIFF_sqrt--STRENGTH 8348 0.4973 < 0.001 0.6463 Errors
## 5 DIFF_log--STRENGTH_sqrt 8380 0.5184 < 0.001 0.6746 Errors
## 6 DIFF_cbrt--STRENGTH 8461 0.0778 0.0251 RT
## 7 DIFF_sqrt--STRENGTH 8480 0.0766 < 0.001 0.0246 RT
## 8 DIFF_cbrt--STRENGTH_sqrt 8520 0.0740 < 0.001 0.0268 RT
## 9 DIFF_log--STRENGTH 8534 0.0734 < 0.001 0.0233 RT
## 10 DIFF_sqrt--STRENGTH_sqrt 8540 0.0728 < 0.001 0.0263 RT
plot_descriptive_err(data, side = "updown")
formula <- brms::bf(
Error ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_poggendorff_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 3086.8 seconds.
## Chain 3 finished in 3150.5 seconds.
## Chain 2 finished in 3172.7 seconds.
## Chain 1 finished in 3187.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 3149.4 seconds.
## Total execution time: 3187.6 seconds.
# parameters::parameters(model_poggendorff_err)
save(model_poggendorff_err, file="models/model_poggendorff_err.Rdata")
load("models/model_poggendorff_err.RData")
p_poggendorff_err <- plot_model_err(data, model_poggendorff_err)
p_poggendorff_err
data <- filter(df, Illusion_Type == "Poggendorff", Error == 0)
plot_descriptive_rt(data, side = "updown")
formula <- brms::bf(
RT ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
# sigma ~ 1 + (1 | Participant),
family = "exgaussian"
)
model_poggendorff_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_ponzo_rt, file="models/model_poggendorff_rt.Rdata")
load("models/model_poggendorff_rt.RData")
plot_ppcheck(model_poggendorff_rt)
# performance::check_model(model_poggendorff_rt)
# parameters::parameters(model)
p_poggendorff_rt <- plot_model_rt(data, model_poggendorff_rt)
p_poggendorff_rt
data <- filter(df, Illusion_Type == "Contrast")
best_models(data)
## Name BIC R2_marginal BF Side Model
## 1 DIFF_sqrt--STRENGTH_sqrt 7837 0.610 4.70e-19 Errors
## 2 DIFF--STRENGTH_sqrt 7842 0.621 = 0.084 4.55e-19 Errors
## 3 DIFF_cbrt--STRENGTH_sqrt 7847 0.607 = 0.008 5.12e-19 Errors
## 4 DIFF_sqrt--STRENGTH_log 7859 0.658 < 0.001 8.41e-19 Errors
## 5 DIFF--STRENGTH_log 7864 0.675 < 0.001 8.08e-19 Errors
## 6 DIFF_sqrt--STRENGTH_sqrt 10454 0.122 2.42e-02 RT
## 7 DIFF_cbrt--STRENGTH_sqrt 10454 0.122 = 0.824 2.43e-02 RT
## 8 DIFF_log--STRENGTH_sqrt 10457 0.122 = 0.211 2.45e-02 RT
## 9 DIFF--STRENGTH_sqrt 10462 0.122 = 0.013 2.44e-02 RT
## 10 DIFF_sqrt--STRENGTH 10479 0.121 < 0.001 1.94e-02 RT
data <- filter(df, Illusion_Type == "Contrast")
plot_descriptive_err(data, side = "updown")
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
family = "bernoulli"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
model_contrast_err <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 3483.4 seconds.
## Chain 2 finished in 3506.1 seconds.
## Chain 4 finished in 3551.5 seconds.
## Chain 3 finished in 3555.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 3524.0 seconds.
## Total execution time: 3555.6 seconds.
# parameters::parameters(model_contrast_err)
save(model_contrast_err, file="models/model_contrast_err.Rdata")
load("models/model_contrast_err.RData")
p_contrast_err <- plot_model_err(data, model_contrast_err)
p_contrast_err
data <- filter(df, Illusion_Type == "Contrast", Error == 0)
plot_descriptive_rt(data, side = "updown")
formula <- brms::bf(
RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
(1 | Participant),
# sigma ~ 1 + (1 | Participant),
family = "exgaussian"
)
model_contrast_rt <- brms::brm(formula,
data = data,
refresh = 0
)
save(model_contrast_rt, file="models/model_contrast_rt.Rdata")
load("models/model_contrast_rt.RData")
plot_ppcheck(model_contrast_rt)
# performance::check_model(model_contrast_rt)
# parameters::parameters(model)
p_contrast_rt <- plot_model_rt(data, model_contrast_rt)
p_contrast_rt
random <- rbind(
extract_random(model_delboeuf_err, "Delboeuf"),
# extract_random(model_delboeuf_rt, "Delboeuf")
extract_random(model_ebbinghaus_err, "Ebbinghaus"),
# extract_random(model_ebbinghaus_rt, "Ebbinghaus"),
extract_random(model_rodframe_err, "Rod-Frame"),
# extract_random(model_rodframe_rt, "Rod-Frame"),
extract_random(model_verticalhorizontal_err, "Vertical-Horizontal"),
# extract_random(model_verticalhorizontal_rt, "Vertical-Horizontal"),
extract_random(model_zollner_err, "Zöllner"),
# extract_random(model_zollner_rt, "Zöllner"),
extract_random(model_white_err, "White"),
# extract_random(model_white_rt, "White"),
extract_random(model_mullerlyer_err, "Müller-Lyer"),
# extract_random(model_mullerlyer_rt, "Müller-Lyer"),
extract_random(model_ponzo_err, "Ponzo"),
# extract_random(model_ponzo_rt, "Ponzo"),
extract_random(model_poggendorff_err, "Poggendorff"),
# extract_random(model_poggendorff_rt, "Poggendorff"),
extract_random(model_contrast_err, "Contrast")
# extract_random(model_contrast_rt, "Contrast")
) |>
filter(
# !str_detect(Parameter, "Intercept_"),
!str_detect(Parameter, "Cong_"),
!str_detect(Parameter, "Null_"))
random |>
# group_by(Illusion_Type, Parameter) |>
# standardize(select="Value") |>
# ungroup() |>
# mutate(order = mean(Value)) |>
# ungroup() |>
# mutate(Participant = fct_reorder(Participant, order)) |>
ggplot(aes(y = Participant, x = Value)) +
# ggdist::stat_slabinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
ggdist::stat_pointinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
# coord_cartesian(xlim=c(-4, 4)) +
facet_wrap(~Parameter, scales = "free_x") +
theme(axis.text.y = element_blank())
random <- random |>
mutate(Parameter = paste0(Illusion_Type, "_", Parameter),
Parameter = format_illusionName(Parameter)) |>
group_by(Participant, Parameter) |>
summarize(Score = median(Value), .groups = "drop") |>
ungroup() |>
group_by(Parameter) |>
mutate(Score = Score / sd(Score)) |> # Scale
ungroup() |>
pivot_wider(names_from = Parameter, values_from = Score)
# plot(estimate_density(select(random, -Participant)))
# Average when incongruent
empirical <- df |>
filter(Illusion_Effect == "Incongruent") |>
group_by(Participant, Illusion_Type) |>
summarize(
# n = n(),
Error_General = sum(Error, na.rm=TRUE) / n(),
RT_Mean_General = mean(RT, na.rm=TRUE),
RT_SD_General = sd(RT, na.rm=TRUE),
) |>
ungroup() |>
mutate(Illusion_Type = format_illusionName(Illusion_Type)) |>
pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_General", "RT_Mean_General", "RT_SD_General")) |>
arrange(Participant)
# Average for 3 higher level of illusion strength
empirical <- df |>
filter(Illusion_Effect == "Incongruent") |>
group_by(Participant, Illusion_Type, Illusion_Strength) |>
summarize(
# n = n(),
Error_Strong = sum(Error, na.rm=TRUE) / n(),
RT_Mean_Strong = mean(RT, na.rm=TRUE),
RT_SD_Strong = sd(RT, na.rm=TRUE),
) |>
slice((n()-2):n()) |>
group_by(Participant, Illusion_Type) |>
summarize(
# n = n(),
Error_Strong = mean(Error_Strong, na.rm=TRUE),
RT_Mean_Strong = mean(RT_Mean_Strong, na.rm=TRUE),
RT_SD_Strong = mean(RT_SD_Strong, na.rm=TRUE),
) |>
ungroup() |>
mutate(Illusion_Type = format_illusionName(Illusion_Type)) |>
pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_Strong", "RT_Mean_Strong", "RT_SD_Strong")) |>
arrange(Participant) |>
cbind(empirical)
# Correlation
r <- data.frame()
for(ill in c("Delboeuf", "Ebbinghaus",
# "RodFrame",
"VerticalHorizontal",
# "Zollner",
# "White",
"MullerLyer",
"Ponzo",
# "Poggendorff",
"Contrast")) {
r <- correlation::correlation(
data = random |>
arrange(Participant) |>
data_select(select=ill, regex=TRUE),
data2 = empirical |>
data_select(select=ill, regex=TRUE)
) |>
mutate(Illusion_Type = ill,
Parameter1 = clean_name(Parameter1),
Parameter2 = clean_name(Parameter2)) |>
select(Illusion_Type, Parameter1, Parameter2, r, p) |>
rbind(r)
}
r |>
ggplot(aes(x = Parameter2, y = Parameter1)) +
geom_tile(aes(fill = r, alpha = p)) +
scale_alpha_continuous(range = c(1, 0)) +
scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
facet_wrap(~Illusion_Type, scales = "free")
cor <- correlation::correlation(random, redundant = TRUE)
p_data <- cor |>
mutate(
Parameter1 = clean_name(Parameter1),
Parameter2 = clean_name(Parameter2)
) |>
correlation::cor_sort(hclust_method = "ward.D2") |>
cor_lower() |>
mutate(
Text = insight::format_value(r, zap_small = TRUE, digits = 3),
Text = str_replace(str_remove(Text, "^0+"), "^-0+", "-"),
Parameter2 = fct_rev(Parameter2)
)
p_data |>
ggplot(aes(x = Parameter2, y = Parameter1)) +
geom_tile(aes(fill = r)) +
# geom_text(aes(label = Text), size = 2) +
scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
labs(title = "Correlation Matrix", x = NULL, y = NULL) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
data <- select(random, -Participant)
rez <- parameters::n_factors(data)
plot(rez)
fa <- parameters::factor_analysis(data,
cor = as.matrix(cor),
n = 8,
rotation = "oblimin",
fm = "mle",
sort = TRUE
)
# psych::omega(data, nfactors=3)
# fa <- psych::fa.multi(as.matrix(cor), nfactors = 3, nfact2 = 1, n.obs = nrow(random))
# psych::fa.multi.diagram(fa)
data <- predict(fa)
rez <- parameters::n_factors(data, rotation = "oblimin")
plot(rez)
fa2 <- parameters::factor_analysis(data,
n = 2,
rotation = "varimax",
fm = "mle",
sort = TRUE
)
# psych::schmid(attributes(fa)$model, nfactors = 1)
library(ggraph)
dat <- rbind(
attributes(fa)$loadings_long,
attributes(fa2)$loadings_long |>
mutate(Component = str_replace(Component, "ML", "G"))
) |>
tidygraph::as_tbl_graph() |>
tidygraph::activate("nodes")
# Get order of indices
# idx <- correlation::cor_sort(cor, hclust_method = "ward.D2") |>
# summary() |>
# pull(Parameter)
# idx <- sort(fa$Variable)
idx <- fa$Variable
get_x <- function(name, x) {
case_when(
name == "G1" ~ 0/1,
name == "G2" ~ 1/1,
name == "ML1" ~ 0/4,
name == "ML2" ~ 1/4,
name == "ML3" ~ 2/4,
name == "ML4" ~ 3/4,
name == "ML5" ~ 4/4,
TRUE ~ x)
}
dat |>
mutate(
x = sapply(as.data.frame(dat)$name, function(x) {
if (x %in% idx) {
(which(idx == x) - 1) / (length(idx) - 1)
} else {
NA
}
})
) |>
mutate(
x = get_x(name, x),
y = case_when(
str_detect(name, "G") ~ 2,
str_detect(name, "ML") ~ 1,
TRUE ~ 0),
angle = ifelse(name %in% idx, 90, 0),
hjust = ifelse(name %in% idx, 1.01, 0),
name = clean_name(name),
Illusion_Type = str_split_fixed(name, " - ", n=2)[, 1]) |>
tidygraph::activate("edges") |>
filter(abs(Loading) > 0.2) |>
ggraph(layout = "nicely") +
geom_edge_link(aes(edge_width = Loading, edge_colour = as.factor(sign(Loading)))) +
geom_node_point(aes(color = Illusion_Type), size=3) +
geom_node_text(aes(label = name, angle = angle, hjust=hjust)) +
scale_y_continuous(expand = expansion(c(.6, .1))) +
scale_edge_width_continuous(range=c(0.1, 1)) +
scale_edge_color_manual(values=c("#F44336", "#2196F3")) +
scale_color_manual(values=c("Delboeuf" = "#F44336", Ebbinghaus = "#E91E63",
"Rod-Frame" = "#9C27B0")) +
guides(edge_colour = "none") +
ggraph::theme_graph()
library(lavaan)
library(tidySEM)
names <- names(select(random, -Participant))
data <- setNames(as.data.frame(str_split_fixed(names, "_", 3)), c("Illusion_Type", "Parameter", "Model"))
data$Name <- names
# Model 1
lvl1 <- paste(data$Illusion_Type, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Illusion_Type))
resid <- as.data.frame(t(utils::combn(x=unique(data$Illusion_Type), 2)))
resid <- paste(resid$V1, "~~", resid$V2)[1:2]
m1 <- sem(model = paste0(c(lvl1, lvl2, resid), collapse = "\n"),
data = random)
# Model 2
lvl1 <- paste(data$Model, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Model))
m2 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
data = random)
# Model 3
lvl1 <- paste(data$Parameter, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Parameter))
m3 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
data = random)
anova(m1, m2, m3)
summary(m1, standardize = TRUE)
tidySEM::graph_sem(model = m1)
# as_tbl_graph(fit, standardize = TRUE)
# Save
dfsub <- dfsub |>
merge(random, by = "Participant") |>
cbind(data, predict(fa2, names=c("G1", "G2")))
write.csv("../data/study3.csv")